function estimation = mesh_estimation(liquid_level, mode_parameters, mode)

vertices = mode_parameters.rotated_vertices;
triangles = mode_parameters.triangles;

%edge_length = 3;%min_bottle_poly_distance(bottle_shape)
%[mode_parameters.vertices mode_parameters.triangles] = triangulation(bottle_shape, edge_length);


%triangle_mass_center = (vertices(triangles(:,1),:) + vertices(triangles(:,2),:) + vertices(triangles(:,3),:)) / 3;
triangles_in_liquid = all([vertices(triangles(:,1),2) <= liquid_level...
                           vertices(triangles(:,2),2) <= liquid_level... 
                           vertices(triangles(:,3),2) <= liquid_level]);
intersected_triangles_index = intersected_triangles_indeces(triangles, vertices, liquid_level);
vertices_number = size(vertices, 1);
triangles_number = size(triangles, 1);
new_vertices = vertices;
new_triangles = triangles;
for i = intersected_triangles_index
    above = [vertices(triangles(i,1),2) > liquid_level...
             vertices(triangles(i,2),2) > liquid_level... 
             vertices(triangles(i,3),2) > liquid_level];
    below = [vertices(triangles(i,1),2) < liquid_level...
             vertices(triangles(i,2),2) < liquid_level... 
             vertices(triangles(i,3),2) < liquid_level];         
    aboveQ = sum(sum(above));
    belowQ = sum(sum(below));
    if (aboveQ == 1 & belowQ == 1)
        t = (liquid_level - vertices(triangles(i, below),2)) /...
            (vertices(triangles(i, below),2) - vertices(triangles(i, above),2)); 
        new_point = vertices(triangles(i, below),:) -...
                   (vertices(triangles(i, above),:) - vertices(triangles(i, below),:)) * t;
        vertices_number = vertices_number + 1;
        new_vertices(vertices_number, :) = new_point;
        new_triangles(i,:) = [vertices_number triangles(i, above) triangles(i, ~above & ~below)];
        triangles_number = triangles_number + 1;
        new_triangles(triangles_number, :) = [vertices_number triangles(i, below) triangles(i, ~above & ~below)];
    else
        if (aboveQ == 2)
            t = (liquid_level - vertices(triangles(i, below),2))./...
                (vertices(triangles(i, below),2) - vertices(triangles(i, above),2));
            mat = vertices(triangles(i, above),:);
            new_point1 = vertices(triangles(i, below),:) -...
                        (mat(1,:) - vertices(triangles(i, below),:)) * t(1);
            new_point2 = vertices(triangles(i, below),:) -...
                        (mat(2,:) - vertices(triangles(i, below),:)) * t(2);        
            new_point3 = (mat(1,:) + mat(2,:)) / 2;        
            vertices_number = vertices_number + 1;
            new_vertices(vertices_number, :) = new_point1;
            vertices_number = vertices_number + 1;
            new_vertices(vertices_number, :) = new_point2;
            vertices_number = vertices_number + 1;
            new_vertices(vertices_number, :) = new_point3;           
            new_triangles(i,:) = [vertices_number - 2 triangles(i, below) vertices_number - 1];
            triangles_number = triangles_number + 1;
            new_triangles(triangles_number, :) = [vertices_number - 2 vertices_number - 1 vertices_number];
            triangles_number = triangles_number + 1;
            above_vert = triangles(i, above);
            new_triangles(triangles_number, :) = [above_vert(1) vertices_number - 2 vertices_number];
            triangles_number = triangles_number + 1;
            new_triangles(triangles_number, :) = [above_vert(2) vertices_number - 1 vertices_number];
        else
            t = (liquid_level - vertices(triangles(i, below),2))./...
                (vertices(triangles(i, below),2) - vertices(triangles(i, above),2));
            mat = vertices(triangles(i, below),:);
            new_point1 = mat(1,:) -...
                        (vertices(triangles(i, above),:) - mat(1,:)) * t(1);
            new_point2 = mat(2,:) -...
                        (vertices(triangles(i, above),:) - mat(2,:)) * t(2);        
            new_point3 = (mat(1,:) + mat(2,:)) / 2;        
            vertices_number = vertices_number + 1;
            new_vertices(vertices_number, :) = new_point1;
            vertices_number = vertices_number + 1;
            new_vertices(vertices_number, :) = new_point2;
            vertices_number = vertices_number + 1;
            new_vertices(vertices_number, :) = new_point3;           
            new_triangles(i,:) = [vertices_number - 2 triangles(i, above) vertices_number - 1];
            triangles_number = triangles_number + 1;
            new_triangles(triangles_number, :) = [vertices_number - 2 vertices_number - 1 vertices_number];
            triangles_number = triangles_number + 1;
            below_vert = triangles(i, below);
            new_triangles(triangles_number, :) = [below_vert(1) vertices_number - 2 vertices_number];
            triangles_number = triangles_number + 1;
            new_triangles(triangles_number, :) = [below_vert(2) vertices_number - 1 vertices_number];
        end
    end
end

% figure;
% for i = 1:size(triangles, 1)
%     plot(vertices(triangles(i,[1 2 3 1]), 1), vertices(triangles(i,[1 2 3 1]), 2), 'k-');
%     hold on
% end
% axis equal


triangle_mass_center = (new_vertices(new_triangles(:,1),:) +...
                        new_vertices(new_triangles(:,2),:) +...
                        new_vertices(new_triangles(:,3),:)) / 3;
new_triangles_in_liquid = triangle_mass_center(:,2) < liquid_level; 
% idx = 1:size(new_triangles,1);
% indeces = idx(new_triangles_in_liquid);
% 
% for i = indeces
%     plot(new_vertices(new_triangles(i,[1 2 3 1]), 1), new_vertices(new_triangles(i,[1 2 3 1]), 2), 'b-');
%     %hold on
% end
% axis equal

edgevec1 = new_vertices(new_triangles(new_triangles_in_liquid,2),:) -...
           new_vertices(new_triangles(new_triangles_in_liquid,1),:);
edgevec2 = new_vertices(new_triangles(new_triangles_in_liquid,3),:) -...
           new_vertices(new_triangles(new_triangles_in_liquid,1),:);        
triangle_volume = cross([edgevec1 zeros(size(edgevec1,1),1)],[edgevec2 zeros(size(edgevec2,1),1)]);
triangle_volume = abs(triangle_volume) / 2;
volume = sum(triangle_volume(:,3));

switch mode
    case 'volume'
        estimation = volume;
    case 'mass_center'
        %triangles_number = sum(sum(new_triangles_in_liquid));
        estimation = (sum([triangle_volume(:,3) triangle_volume(:,3)].* triangle_mass_center(new_triangles_in_liquid,:)) / volume)'; 
        %plot(estimation(1), estimation(2), 'bx', 'LineWidth',2, 'MarkerSize', 10);
    otherwise 
        error('Invalid mode');
end



    



end
